/*
 * topological-hex, a program to compute combinatorial hexahedral meshes.
 *
 * Copyright (C) <2018> <Université catholique de Louvain (UCL), Belgique>
 *
 * List of the contributors to the development, description and complete
 * License: see LICENSE file.
 *
 * This program (topological-hex) is free software:
 * you can redistribute it and/or modify it under the terms
 * of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program (see COPYING file).  If not,
 * see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "solver.h"
#include "shellable_dfs.h"
#include "symmetry.h"
#include "lower_bound.h"
#include "dual.h"

static
quad_edge *other_side(solver *solver, quad_edge e) {
  return &solver->quad_adjacency[4 * quad_edge_quad(e) + quad_edge_edge(e)];
}

static
vertex_index start_v(const solver *solver, quad_edge e) {
  return solver->quads[4 * quad_edge_quad(e) + quad_edge_edge(e)];
}

static
vertex_index end_v(const solver *solver, quad_edge e) {
  return solver->quads[4 * quad_edge_quad(e) + ((quad_edge_edge(e) + 1) % 4)];
}

static
hex_index hex_containing(const solver *solver, quad_index q) {
  vertex_index a = solver->quads[4 * q + 0];
  vertex_index b = solver->quads[4 * q + 2];

  return quad_descriptor_hex_index(
    solver->quad_descriptors[diagonal_index(solver->num_vertices, a, b)]);
}

static
bool check_edge(const solver *solver, vertex_index v_a, vertex_index v_b) {
  return !vertex_set_contains(&solver->known_diagonals[v_a], v_b) &&
    !vertex_set_contains(&solver->known_neighbors[v_a], v_b);
}


static
bool check_diagonal(const solver *solver, vertex_index v_a, vertex_index v_b) {
  return !vertex_set_contains(&solver->known_diagonals[v_a], v_b) &&
    !vertex_set_contains(&solver->known_neighbors[v_a], v_b);
}

static
bool check_quad_diagonal(const solver *solver, vertex_index v_a, vertex_index v_b) {
  return solver->quad_descriptors[diagonal_index(solver->num_vertices, v_a, v_b)]
    == EMPTY_QUAD_DESCRIPTOR;
}

bool allowed_flip_1_to_5(solver *solver, quad_index q) {
  return solver->max_vertex_id + 4 <= solver->num_vertices &&
    solver->max_quad_id + 5 <= solver->num_quads;
}

bool allowed_flip_5_to_1(solver *solver, quad_index q) {
  if (solver->max_quad_id + 1 > solver->num_quads)
    return false;

  hex_index hexes[4];
  quad_index quads[4];
  vertex_index v_out[4];

  hex_index center_hex = hex_containing(solver, q);

  for (quad_edge_index i = 0; i < 4; i++) {
    quad_edge edge = make_quad_edge(q, i);
    quad_edge other_edge = *other_side(solver, edge);

    hexes[i] = hex_containing(solver, quad_edge_quad(other_edge));
    quads[i] = quad_edge_quad(*other_side(solver, opposite_edge(other_edge)));
    v_out[i] = end_v(solver, next_edge(other_edge));

    if (*other_side(
          solver, next_edge(*other_side(solver, next_edge(other_edge)))) !=
        prev_edge(edge))
      return false;

    if (hexes[i] != NO_HEXAHEDRON && hexes[i] == center_hex)
      return false;

    for (quad_vertex_index j = 0; j < i; j++) {
      if (hexes[i] != NO_HEXAHEDRON && hexes[j] == hexes[i]) return false;
      if (quads[j] == quads[i]) return false;
    }
  }

  const vertex_index *v_in = solver->quads + 4 * q;

  for (quad_vertex_index i = 0; i < 4; i++) {
    if (!check_diagonal(solver, v_in[i], v_out[(i + 2) % 4]))
      return false;

    if (!check_quad_diagonal(solver, v_out[i], v_out[(i + 2) % 4]))
      return false;
  }

  return true;
}

bool allowed_flip_2_to_4(solver *solver, quad_edge edge) {
  if (solver->max_vertex_id + 2 > solver->num_vertices)
    return false;

  if (solver->max_quad_id + 4 > solver->num_quads)
    return false;

  quad_edge other_edge = *other_side(solver, edge);

  quad_index q_a = quad_edge_quad(edge);
  quad_index q_b = quad_edge_quad(other_edge);

  if (q_a > q_b) return false; /* to break symmetry */

  hex_index h_a = hex_containing(solver, q_a);
  hex_index h_b = hex_containing(solver, q_b);

  if (h_a != NO_HEXAHEDRON && h_a == h_b) return false;

  quad_index next_a = quad_edge_quad(*other_side(solver, next_edge(edge)));
  quad_index prev_a = quad_edge_quad(*other_side(solver, prev_edge(edge)));
  quad_index next_b = quad_edge_quad(*other_side(solver, next_edge(other_edge)));
  quad_index prev_b = quad_edge_quad(*other_side(solver, prev_edge(other_edge)));

  if (next_a == prev_b || next_b == prev_a) return false;

  if (!check_quad_diagonal(
        solver,
        start_v(solver, opposite_edge(edge)),
        start_v(solver, prev_edge(other_edge))))
    return false;

  if (!check_quad_diagonal(
        solver,
        start_v(solver, opposite_edge(other_edge)),
        start_v(solver, prev_edge(edge))))
    return false;

  if (!check_diagonal(
        solver,
        start_v(solver, opposite_edge(edge)),
        start_v(solver, opposite_edge(other_edge))))
    return false;

  if (!check_diagonal(
        solver,
        start_v(solver, prev_edge(edge)),
        start_v(solver, prev_edge(other_edge))))
    return false;

  return true;
}

bool allowed_flip_4_to_2(solver *solver, quad_edge edge) {
  if (solver->max_quad_id + 2 > solver->num_quads)
    return false;

  quad_edge other_edge = *other_side(solver, edge);

  quad_index q_a = quad_edge_quad(edge);
  quad_index q_b = quad_edge_quad(other_edge);

  hex_index h_a = hex_containing(solver, q_a);
  hex_index h_b = hex_containing(solver, q_b);

  if (h_a != NO_HEXAHEDRON && h_a == h_b) return false;

  if (q_a > q_b) return false; /* to break symmetry */

  if (*other_side(solver, next_edge(*other_side(solver, next_edge(edge))))
      != prev_edge(other_edge))
    return false;

  if (*other_side(solver, next_edge(*other_side(solver, next_edge(other_edge))))
      != prev_edge(edge))
    return false;

  hex_index h_c = hex_containing(
    solver, quad_edge_quad(*other_side(solver, next_edge(edge))));
  hex_index h_d = hex_containing(
    solver, quad_edge_quad(*other_side(solver, next_edge(other_edge))));

  if (h_c != NO_HEXAHEDRON && (h_a == h_c || h_b == h_c)) return false;
  if (h_d != NO_HEXAHEDRON && (h_a == h_d || h_b == h_d || h_c == h_d)) return false;

  quad_index top_q_a = quad_edge_quad(
    *other_side(solver, prev_edge(*other_side(solver, next_edge(edge)))));
  quad_index top_q_b = quad_edge_quad(
    *other_side(solver, opposite_edge(edge)));
  quad_index top_q_c = quad_edge_quad(
    *other_side(solver, next_edge(*other_side(solver, prev_edge(edge)))));

  if (top_q_a == top_q_b || top_q_b == top_q_c || top_q_c == top_q_a)
    return false;

  quad_index bot_q_a = quad_edge_quad(
    *other_side(solver, prev_edge(*other_side(solver, next_edge(other_edge)))));
  quad_index bot_q_b = quad_edge_quad(
    *other_side(solver, opposite_edge(other_edge)));
  quad_index bot_q_c = quad_edge_quad(
    *other_side(solver, next_edge(*other_side(solver, prev_edge(other_edge)))));

  if (bot_q_a == bot_q_b || bot_q_b == bot_q_c || bot_q_c == bot_q_a)
    return false;

  /*
   *  top A ___________ top B
   *       /\         /\
   *      /  \       /  \
   *   A / n a\_____/ n b\ B
   *     \    /     \    /
   *      \  /       \  /
   * bot A \/_________\/ bot B
   */

  vertex_index v_a = start_v(solver, opposite_edge(*other_side(solver, prev_edge(edge))));
  vertex_index v_b = start_v(solver, opposite_edge(*other_side(solver, prev_edge(other_edge))));

  vertex_index top_a = start_v(solver, prev_edge(edge));
  vertex_index top_b = start_v(solver, opposite_edge(edge));

  vertex_index bot_a = start_v(solver, opposite_edge(other_edge));
  vertex_index bot_b = start_v(solver, prev_edge(other_edge));

  vertex_index n_a = start_v(solver, edge);
  vertex_index n_b = start_v(solver, other_edge);

  if (!(check_edge(solver, v_a, v_b) &&

        check_quad_diagonal(solver, top_a, v_b) &&
        check_quad_diagonal(solver, top_b, v_a) &&
        check_quad_diagonal(solver, bot_a, v_b) &&
        check_quad_diagonal(solver, bot_b, v_a) &&

        check_diagonal(solver, top_a, bot_b) &&
        check_diagonal(solver, top_b, bot_a) &&
        check_diagonal(solver, v_b, n_a) &&
        check_diagonal(solver, v_a, n_b))) {
    return false;
  }

  return true;
}

bool allowed_flip_3_line(solver *solver, quad_edge edge_a) {
  if (solver->max_quad_id + 3 > solver->num_quads)
    return false;

  quad_edge other_a = *other_side(solver, edge_a);
  quad_edge edge_b  = opposite_edge(other_a);
  quad_edge other_b = *other_side(solver, edge_b);

  if (quad_edge_quad(edge_a) > quad_edge_quad(other_b))
    return false;

  hex_index hex_a = hex_containing(solver, quad_edge_quad(edge_a));
  hex_index hex_b = hex_containing(solver, quad_edge_quad(edge_b));
  hex_index hex_c = hex_containing(solver, quad_edge_quad(other_b));

  if (hex_a != NO_HEXAHEDRON && hex_a == hex_c) return false;
  if (hex_b != NO_HEXAHEDRON && hex_a == hex_b) return false;
  if (hex_c != NO_HEXAHEDRON && hex_c == hex_b) return false;

  {
    quad_index q0 = quad_edge_quad(*other_side(solver, next_edge(edge_a)));
    quad_index q1 = quad_edge_quad(*other_side(solver, next_edge(edge_b)));
    quad_index q2 = quad_edge_quad(*other_side(solver, prev_edge(other_b)));

    if (q0 == q1 || q1 == q2 || q2 == q0) return false;
  }

  {
    quad_index q0 = quad_edge_quad(*other_side(solver, opposite_edge(edge_a)));
    quad_index q1 = quad_edge_quad(*other_side(solver, opposite_edge(other_b)));

    if (q0 == q1) return false;
  }

  {
    quad_index q0 = quad_edge_quad(*other_side(solver, prev_edge(edge_a)));
    quad_index q1 = quad_edge_quad(*other_side(solver, prev_edge(edge_b)));
    quad_index q2 = quad_edge_quad(*other_side(solver, next_edge(other_b)));

    if (q0 == q1 || q1 == q2 || q2 == q0) return false;
  }

  if (!check_edge(solver,
                  start_v(solver, opposite_edge(other_b)),
                  start_v(solver, prev_edge(edge_a)))) {
    return false;
  }

  if (!check_edge(solver,
                  start_v(solver, opposite_edge(edge_a)),
                  start_v(solver, prev_edge(other_b)))) {
    return false;
  }

  if (!check_diagonal(solver,
                      start_v(solver, opposite_edge(edge_a)),
                      start_v(solver, next_edge(other_b)))) {
    return false;
  }

  if (!check_diagonal(solver,
                      start_v(solver, next_edge(edge_a)),
                      start_v(solver, opposite_edge(other_b)))) {
    return false;
  }

  if (!check_diagonal(solver,
                      start_v(solver, other_b),
                      start_v(solver, prev_edge(edge_a)))) {
    return false;
  }

  if (!check_diagonal(solver,
                      start_v(solver, edge_a),
                      start_v(solver, prev_edge(other_b)))) {
    return false;
  }

  if (!check_quad_diagonal(solver,
                           start_v(solver, edge_a),
                           start_v(solver, opposite_edge(other_b)))) {
    return false;
  }

  if (!check_quad_diagonal(solver,
                           start_v(solver, prev_edge(edge_a)),
                           start_v(solver, next_edge(other_b)))) {
    return false;
  }

  if (!check_quad_diagonal(solver,
                           start_v(solver, next_edge(edge_a)),
                           start_v(solver, prev_edge(other_b)))) {
    return false;
  }

  if (!check_quad_diagonal(solver,
                           start_v(solver, other_b),
                           start_v(solver, opposite_edge(edge_a)))) {
    return false;
  }

  if (!check_quad_diagonal(solver,
                           start_v(solver, prev_edge(edge_a)),
                           start_v(solver, prev_edge(other_b)))) {
    return false;
  }

  if (!check_quad_diagonal(solver,
                           start_v(solver, opposite_edge(edge_a)),
                           start_v(solver, opposite_edge(other_b)))) {
    return false;
  }

  return true;
}

bool allowed_flip_3_cycle(solver *solver, quad_edge edge_a) {
  if (solver->max_vertex_id + 1 > solver->num_vertices)
    return false;

  if (solver->max_quad_id + 3 > solver->num_quads)
    return false;

  quad_edge other_a = *other_side(solver, edge_a);

  quad_edge edge_b = next_edge(other_a);
  quad_edge other_b = *other_side(solver, edge_b);

  quad_edge edge_c = next_edge(other_b);
  quad_edge other_c = *other_side(solver, edge_c);

  quad_index q_a = quad_edge_quad(edge_a);
  quad_index q_b = quad_edge_quad(edge_b);
  quad_index q_c = quad_edge_quad(edge_c);

  if (quad_edge_quad(other_c) != q_a)
    return false;

  if (q_a >= q_c || q_a >= q_b || q_b == q_c)
    return false;

  hex_index hex_a = hex_containing(solver, q_a);
  hex_index hex_b = hex_containing(solver, q_b);
  hex_index hex_c = hex_containing(solver, q_c);

  if (hex_a != NO_HEXAHEDRON && (hex_a == hex_b || hex_a == hex_c))
    return false;

  if (hex_b != NO_HEXAHEDRON && hex_b == hex_c)
    return false;

  if (quad_edge_quad(*other_side(solver, next_edge(edge_a))) ==
      quad_edge_quad(*other_side(solver, prev_edge(other_a)))) {
    return false;
  }

  if (quad_edge_quad(*other_side(solver, next_edge(edge_b))) ==
      quad_edge_quad(*other_side(solver, prev_edge(other_b)))) {
    return false;
  }

  if (quad_edge_quad(*other_side(solver, next_edge(edge_c))) ==
      quad_edge_quad(*other_side(solver, prev_edge(other_c)))) {
    return false;
  }

  if (!check_diagonal(solver,
                      start_v(solver, opposite_edge(edge_a)),
                      start_v(solver, next_edge(edge_b)))) {
    return false;
  }

  if (!check_diagonal(solver,
                      start_v(solver, opposite_edge(edge_b)),
                      start_v(solver, next_edge(edge_c)))) {
    return false;
  }

  if (!check_diagonal(solver,
                      start_v(solver, opposite_edge(edge_c)),
                      start_v(solver, next_edge(edge_a)))) {
    return false;
  }

  if (!check_quad_diagonal(solver,
                           start_v(solver, opposite_edge(edge_a)),
                           start_v(solver, opposite_edge(edge_b)))) {
    return false;
  }

  if (!check_quad_diagonal(solver,
                           start_v(solver, opposite_edge(edge_b)),
                           start_v(solver, opposite_edge(edge_c)))) {
    return false;
  }

  if (!check_quad_diagonal(solver,
                           start_v(solver, opposite_edge(edge_c)),
                           start_v(solver, opposite_edge(edge_a)))) {
    return false;
  }

  return true;
}

bool allowed_finish_mesh(solver *solver, quad_index *boundary) {
  quad_index num_quads = 0;
  hex_index owners[6];

  for (hex_index i = solver->quad_queue_offset;
       i < solver->quad_queue_size; i++) {
    if (solver->num_occurrences[solver->quad_queue[i]] != 0) {
      if (num_quads == 6) return false;
      boundary[num_quads] = solver->quad_queue[i];
      owners[num_quads]   = hex_containing(solver, solver->quad_queue[i]);
      num_quads++;
    }
  }

  if (num_quads != 6)
    return false;

  for (hex_facet_index i = 0; i < 6; i++) {
    if (owners[i] != NO_HEXAHEDRON) {
      for (hex_facet_index j = 0; j < i; j++) {
        if (owners[j] == owners[i]) return false;
      }
    }
  }

  vertex_index hex[8];
  memcpy(hex, solver->quads + 4 * boundary[0], 4 * sizeof(*hex));

  for (quad_edge_index i = 0; i < 4; i++) {
    quad_edge edge = make_quad_edge(boundary[0], i);
    quad_edge other_edge = *other_side(solver, edge);
    hex[4 + i] = start_v(solver, opposite_edge(other_edge));
  }

  return check_diagonal(solver, hex[0], hex[6]) &&
    check_diagonal(solver, hex[1], hex[7]) &&
    check_diagonal(solver, hex[2], hex[4]) &&
    check_diagonal(solver, hex[3], hex[5]);
}

void finish_mesh(solver *solver, const quad_index *boundary) {
  vertex_index hex[8];
  memcpy(hex, solver->quads + 4 * boundary[0], 4 * sizeof(*hex));

  quad_index bottom = boundary[0];

  for (quad_edge_index i = 0; i < 4; i++) {
    quad_edge edge = make_quad_edge(bottom, i);
    quad_edge other_edge = *other_side(solver, edge);
    hex[4 + i] = start_v(solver, opposite_edge(other_edge));
  }

  /* for (int i = 0; i < 2; i++) { */
  /*   vertex_index t = hex[4*i + 1]; */
  /*   hex[4*i + 1] = hex[4*i + 3]; */
  /*   hex[4*i + 3] = t; */
  /* } */

  memcpy(solver->hexes + 8 * solver->hex_id, hex, sizeof(hex));

  for (hex_facet_index i = 0; i < 6; i++)
    solver->num_occurrences[boundary[i]]--;
}

void undo_finish_mesh(solver *solver, const quad_index *boundary) {
  for (hex_facet_index i = 0; i < 6; i++)
    solver->num_occurrences[boundary[i]]++;
}

static
void create_shellable_delta(solver *solver, shellable_delta *delta,
                            uint8_t num_new_vertices,
                            uint8_t num_finished_vertices, const vertex_index *vertices,
                            uint8_t num_new_quads, const quad_index *quads) {
  delta->num_new_vertices = num_new_vertices;

  solver->max_vertex_id += num_new_vertices;

  uint8_t num_new_edges = 0;
  for (hex_vertex_index i = 0; i < 8; i++) {
    vertex_index a = solver->hexes[8 * solver->hex_id + i];

    for (uint8_t j = 0; j < 3; j++) {
      if (HEX_NEIGHBORS[i][j] < i) {
        vertex_index b = solver->hexes[8 * solver->hex_id + HEX_NEIGHBORS[i][j]];
        if (!vertex_set_contains(&solver->known_neighbors[a], b)) {
          vertex_set_insert(&solver->known_neighbors[a], b);
          vertex_set_insert(&solver->known_neighbors[b], a);

          delta->new_edges[2 * num_new_edges + 0] = a;
          delta->new_edges[2 * num_new_edges + 1] = b;
          num_new_edges++;
        }
      }
    }
  }

  delta->num_new_edges = num_new_edges;

  const vertex_index *hex = solver->hexes + 8 * solver->hex_id;
  for (quad_vertex_index i = 0; i < 4; i++) {
    vertex_index a = hex[i];
    vertex_index b = hex[4 + (i + 2) % 4];
    vertex_set_insert(&solver->known_diagonals[a], b);
    vertex_set_insert(&solver->known_diagonals[b], a);
  }

  delta->num_new_quads = num_new_quads;
  solver->max_quad_id += num_new_quads;

  uint8_t num_old_quads = 6 - num_new_quads;
  memcpy(delta->quads, quads, num_old_quads * sizeof(*quads));

  for (uint8_t i = num_old_quads; i < 6; i++) {
    solver->quad_queue[solver->quad_queue_size++] = quads[i];

    const vertex_index *vs = solver->quads + 4 * quads[i];
    solver->quad_descriptors[diagonal_index(solver->num_vertices, vs[0], vs[2])] =
      make_quad_descriptor(quads[i], solver->hex_id, vs[1], vs[3]);
    solver->quad_descriptors[diagonal_index(solver->num_vertices, vs[1], vs[3])] =
      make_quad_descriptor(quads[i], solver->hex_id, vs[0], vs[2]);

    for (quad_edge_index j = 0; j < 4; j++) {
      *other_side(solver, solver->quad_adjacency[4 * quads[i] + j]) =
        make_quad_edge(quads[i], j);
    }
  }

  for (hex_facet_index i = 0; i < 6; i++)
    solver->num_occurrences[quads[i]]--;

  delta->queue_offset_change = 0;
  while (solver->quad_queue_offset < solver->quad_queue_size &&
         solver->num_occurrences[solver->quad_queue[solver->quad_queue_offset]] == 0) {
    delta->queue_offset_change++;
    solver->quad_queue_offset++;
  }

  delta->num_finished_vertices = num_finished_vertices;
  for (uint8_t i = 0; i < num_finished_vertices; i++) {
    delta->finished_vertices[i] = vertices[i];
    vertex_set_insert(&solver->finished_vertices, vertices[i]);
  }
}

void undo_shellable_delta(solver *solver, const shellable_delta *delta) {
  uint8_t num_old_quads = 6 - delta->num_new_quads;

  solver->max_vertex_id -= delta->num_new_vertices;

  for (uint8_t i = 0; i < delta->num_new_quads; i++) {
    quad_index q = solver->quad_queue[solver->quad_queue_size - 1];
    solver->num_occurrences[q]++;

    const vertex_index *vs = solver->quads + 4 * q;
    solver->quad_descriptors[diagonal_index(solver->num_vertices, vs[0], vs[2])] =
      EMPTY_QUAD_DESCRIPTOR;
    solver->quad_descriptors[diagonal_index(solver->num_vertices, vs[1], vs[3])] =
      EMPTY_QUAD_DESCRIPTOR;

    solver->quad_queue_size--;
  }

  solver->max_quad_id -= delta->num_new_quads;

  for (uint8_t i = 0; i < num_old_quads; i++) {
    quad_index q = delta->quads[i];
    solver->num_occurrences[q]++;
    for (quad_edge_index j = 0; j < 4; j++) {
      *other_side(solver, solver->quad_adjacency[4 * q + j]) =
        make_quad_edge(q, j);
    }
  }

  solver->quad_queue_offset -= delta->queue_offset_change;

  for (uint8_t i = 0; i < delta->num_new_edges; i++) {
    vertex_set_remove(&solver->known_neighbors[delta->new_edges[2 * i + 0]],
                      delta->new_edges[2 * i + 1]);
    vertex_set_remove(&solver->known_neighbors[delta->new_edges[2 * i + 1]],
                      delta->new_edges[2 * i + 0]);
  }

  const vertex_index *hex = solver->hexes + 8 * solver->hex_id;
  for (quad_vertex_index i = 0; i < 4; i++) {
    vertex_index a = hex[i];
    vertex_index b = hex[4 + (i + 2) % 4];
    vertex_set_remove(&solver->known_diagonals[a], b);
    vertex_set_remove(&solver->known_diagonals[b], a);
  }

  for (uint8_t i = 0; i < delta->num_finished_vertices; i++) {
    vertex_set_remove(&solver->finished_vertices,
                      delta->finished_vertices[i]);
  }
}

static
size_t recurse_shellable_dfs(solver *solver, shellable_delta *delta,
                             no_goods_delta *symmetry_delta,
                             const uint8_t *trail) {
  solver->bound_mask = ~(uint32_t)0;

  if (incremental_domination_test(solver)) {
    undo_shellable_delta(solver, delta);
    incremental_symmetry_data_restore(solver, trail);
    return 0;
  }
  else {
    shellable_dfs(solver);
    undo_shellable_delta(solver, delta);
    incremental_symmetry_data_restore(solver, trail);
    solver_hex_is_no_good(solver, symmetry_delta);
    return 1;
  }
}

void flip_1_to_5(solver *solver, quad_index quad, shellable_delta *delta) {
  const vertex_index *bot = solver->quads + 4 * quad;
  const vertex_index top[4] = {
    solver->max_vertex_id + 0,
    solver->max_vertex_id + 1,
    solver->max_vertex_id + 2,
    solver->max_vertex_id + 3,
  };

  memcpy(solver->hexes + 8 * solver->hex_id, bot, 4 * sizeof(*bot));
  memcpy(solver->hexes + 8 * solver->hex_id + 4, top, 4 * sizeof(*top));

  quad_index quads[] = {
    quad,

    solver->max_quad_id + 0,
    solver->max_quad_id + 1,
    solver->max_quad_id + 2,
    solver->max_quad_id + 3,
    solver->max_quad_id + 4,
  };

  const quad_edge *neighbors = solver->quad_adjacency + 4 * quad;

  for (quad_edge_index i = 0; i < 4; i++) {
    vertex_index new_quad[4] = {
      bot[i], bot[(i + 1) % 4], top[(i + 1) % 4], top[i]
    };

    quad_edge adj[4] = {
      neighbors[i],
      make_quad_edge(quads[1 + ((i + 1) % 4)], 3),
      make_quad_edge(quads[5], i),
      make_quad_edge(quads[1 + ((i + 3) % 4)], 1),
    };

    memcpy(solver->quads + 4 * quads[1 + i], new_quad, sizeof(new_quad));
    memcpy(solver->quad_adjacency + 4 * quads[1 + i], adj, sizeof(adj));
  }

  quad_edge top_adj[4] = {
    make_quad_edge(quads[1 + 1], 2),
    make_quad_edge(quads[1 + 2], 2),
    make_quad_edge(quads[1 + 3], 2),
    make_quad_edge(quads[1 + 4], 2),
  };

  memcpy(solver->quads + 4 * quads[5], top, sizeof(top));
  memcpy(solver->quad_adjacency + 4 * quads[5], top_adj, sizeof(top_adj));

  create_shellable_delta(solver, delta, 4, 0, NULL, 5, quads);
}

void flip_5_to_1(solver *solver, quad_index q, shellable_delta *delta) {
  quad_index quads[6] = {
    q,
    quad_edge_quad(*other_side(solver, make_quad_edge(q, 0))),
    quad_edge_quad(*other_side(solver, make_quad_edge(q, 1))),
    quad_edge_quad(*other_side(solver, make_quad_edge(q, 2))),
    quad_edge_quad(*other_side(solver, make_quad_edge(q, 3))),

    solver->max_quad_id,
  };

  vertex_index new_quad[4] = {
    start_v(solver, opposite_edge(*other_side(solver, make_quad_edge(q, 0)))),
    start_v(solver, opposite_edge(*other_side(solver, make_quad_edge(q, 1)))),
    start_v(solver, opposite_edge(*other_side(solver, make_quad_edge(q, 2)))),
    start_v(solver, opposite_edge(*other_side(solver, make_quad_edge(q, 3)))),
  };

  quad_edge new_adj[4] = {
    *other_side(solver, opposite_edge(*other_side(solver, make_quad_edge(q, 0)))),
    *other_side(solver, opposite_edge(*other_side(solver, make_quad_edge(q, 1)))),
    *other_side(solver, opposite_edge(*other_side(solver, make_quad_edge(q, 2)))),
    *other_side(solver, opposite_edge(*other_side(solver, make_quad_edge(q, 3)))),
  };

  memcpy(solver->hexes + 8 * solver->hex_id, solver->quads + 4 * q,  4 * sizeof(*solver->hexes));
  memcpy(solver->hexes + 8 * solver->hex_id + 4, new_quad, sizeof(new_quad));

  memcpy(solver->quads + 4 * solver->max_quad_id, new_quad, sizeof(new_quad));
  memcpy(solver->quad_adjacency + 4 * solver->max_quad_id, new_adj, sizeof(new_adj));

  create_shellable_delta(solver, delta, 0,
                         4, solver->quads + 4 * q,
                         1, quads);
}

void flip_2_to_4(solver *solver, quad_edge edge, shellable_delta *delta) {
  quad_edge other_edge = *other_side(solver, edge);

  vertex_index hex[8] = {
    start_v(solver, edge),
    start_v(solver, next_edge(edge)),
    start_v(solver, opposite_edge(edge)),
    start_v(solver, prev_edge(edge)),

    start_v(solver, opposite_edge(other_edge)),
    start_v(solver, prev_edge(other_edge)),
    solver->max_vertex_id + 0,
    solver->max_vertex_id + 1,
  };

  quad_index quads[6] = {
    quad_edge_quad(edge),
    quad_edge_quad(other_edge),

    solver->max_quad_id + 0,
    solver->max_quad_id + 1,
    solver->max_quad_id + 2,
    solver->max_quad_id + 3,
  };

  vertex_index quad_data[4][4] = {
    {
      solver->max_vertex_id + 1, solver->max_vertex_id + 0,
      start_v(solver, opposite_edge(edge)),
      start_v(solver, prev_edge(edge)),
    },

    {
      solver->max_vertex_id + 0, solver->max_vertex_id + 1,
      start_v(solver, opposite_edge(other_edge)),
      start_v(solver, prev_edge(other_edge)),
    },

    {
      solver->max_vertex_id + 0,
      start_v(solver, prev_edge(other_edge)),
      start_v(solver, next_edge(edge)),
      start_v(solver, opposite_edge(edge)),
    },

    {
      solver->max_vertex_id + 1,
      start_v(solver, prev_edge(edge)),
      start_v(solver, next_edge(other_edge)),
      start_v(solver, opposite_edge(other_edge)),
    }
  };

  quad_edge adj_data[4][4] = {
    {
      make_quad_edge(quads[2 + 1], 0),
      make_quad_edge(quads[2 + 2], 3),
      *other_side(solver, opposite_edge(edge)),
      make_quad_edge(quads[2 + 3], 0),
    },

    {
      make_quad_edge(quads[2 + 0], 0),
      make_quad_edge(quads[2 + 3], 3),
      *other_side(solver, opposite_edge(other_edge)),
      make_quad_edge(quads[2 + 2], 0),
    },

    {
      make_quad_edge(quads[2 + 0], 3),
      *other_side(solver, prev_edge(other_edge)),
      *other_side(solver, next_edge(edge)),
      make_quad_edge(quads[2 + 1], 1),
    },

    {
      make_quad_edge(quads[2 + 0], 3),
      *other_side(solver, prev_edge(edge)),
      *other_side(solver, next_edge(other_edge)),
      make_quad_edge(quads[2 + 1], 1),
    },
  };

  memcpy(solver->hexes + 8 * solver->hex_id, hex, sizeof(hex));
  memcpy(solver->quads + 4 * solver->max_quad_id,
         quad_data, sizeof(quad_data));
  memcpy(solver->quad_adjacency + 4 * solver->max_quad_id,
         adj_data, sizeof(adj_data));

  create_shellable_delta(solver, delta, 2, 0, NULL, 4, quads);
}

void flip_4_to_2(solver *solver, quad_edge edge, shellable_delta *delta) {
  quad_edge other = *other_side(solver, edge);

  quad_index quads[6] = {
    quad_edge_quad(edge),
    quad_edge_quad(other),
    quad_edge_quad(*other_side(solver, next_edge(edge))),
    quad_edge_quad(*other_side(solver, next_edge(other))),

    solver->max_quad_id + 0,
    solver->max_quad_id + 1,
  };

  vertex_index quad_data[2][4] = {
    {
      start_v(solver, prev_edge(*other_side(solver, next_edge(other)))),
      start_v(solver, prev_edge(*other_side(solver, next_edge(edge)))),
      start_v(solver, opposite_edge(edge)),
      start_v(solver, prev_edge(edge)),
    },

    {
      start_v(solver, prev_edge(*other_side(solver, next_edge(edge)))),
      start_v(solver, prev_edge(*other_side(solver, next_edge(other)))),
      start_v(solver, opposite_edge(other)),
      start_v(solver, prev_edge(other)),
    }
  };

  quad_edge adj_data[2][4] = {
    {
      make_quad_edge(quads[4 + 1], 0),
      *other_side(solver, prev_edge(*other_side(solver, next_edge(edge)))),
      *other_side(solver, opposite_edge(edge)),
      *other_side(solver, next_edge(*other_side(solver, prev_edge(edge)))),
    },

    {
      make_quad_edge(quads[4 + 0], 0),
      *other_side(solver, prev_edge(*other_side(solver, next_edge(other)))),
      *other_side(solver, opposite_edge(other)),
      *other_side(solver, next_edge(*other_side(solver, prev_edge(other)))),
    }
  };

  vertex_index hex[8] = {
    quad_data[1][2], start_v(solver, edge), quad_data[0][3], quad_data[0][0],
    quad_data[1][3], start_v(solver, other), quad_data[0][2], quad_data[0][1],
  };

  vertex_index finished[2] = {start_v(solver, edge), end_v(solver, edge)};

  memcpy(solver->hexes + 8 * solver->hex_id, hex, sizeof(hex));
  memcpy(solver->quads + 4 * solver->max_quad_id,
         quad_data, sizeof(quad_data));
  memcpy(solver->quad_adjacency + 4 * solver->max_quad_id,
         adj_data, sizeof(adj_data));

  create_shellable_delta(solver, delta, 0, 2, finished, 2, quads);
}

void flip_3_line(solver *solver, quad_edge edge_a, shellable_delta *delta) {
  quad_edge other_a = *other_side(solver, edge_a);
  quad_edge edge_b  = opposite_edge(other_a);
  quad_edge other_b = *other_side(solver, edge_b);

  quad_index quads[6] = {
    quad_edge_quad(edge_a),
    quad_edge_quad(edge_b),
    quad_edge_quad(other_b),

    solver->max_quad_id + 0,
    solver->max_quad_id + 1,
    solver->max_quad_id + 2,
  };

  vertex_index hex[8] = {
    start_v(solver, edge_a),
    start_v(solver, next_edge(edge_a)),
    start_v(solver, opposite_edge(edge_a)),
    start_v(solver, prev_edge(edge_a)),

    start_v(solver, next_edge(other_b)),
    start_v(solver, other_b),
    start_v(solver, prev_edge(other_b)),
    start_v(solver, opposite_edge(other_b)),
  };

  quad_index quad_data[3][4] = {
    {
      start_v(solver, opposite_edge(edge_a)),
      start_v(solver, prev_edge(other_b)),
      start_v(solver, other_b),
      start_v(solver, next_edge(edge_a)),
    },

    {
      start_v(solver, prev_edge(other_b)),
      start_v(solver, opposite_edge(edge_a)),
      start_v(solver, prev_edge(edge_a)),
      start_v(solver, opposite_edge(other_b)),
    },

    {
      start_v(solver, opposite_edge(other_b)),
      start_v(solver, prev_edge(edge_a)),
      start_v(solver, edge_a),
      start_v(solver, edge_b),
    }
  };

  quad_edge adj_data[3][4] = {
    {
      make_quad_edge(quads[3 + 1], 0),
      *other_side(solver, prev_edge(other_b)),
      *other_side(solver, prev_edge(other_a)),
      *other_side(solver, next_edge(edge_a)),
    },

    {
      make_quad_edge(quads[3 + 1], 0),
      *other_side(solver, opposite_edge(edge_a)),
      make_quad_edge(quads[3 + 2], 0),
      *other_side(solver, opposite_edge(other_b)),
    },

    {
      make_quad_edge(quads[3 + 1], 2),
      *other_side(solver, prev_edge(edge_a)),
      *other_side(solver, next_edge(other_a)),
      *other_side(solver, next_edge(other_b)),
    }
  };

  memcpy(solver->hexes + 8 * solver->hex_id, hex, sizeof(hex));
  memcpy(solver->quads + 4 * solver->max_quad_id,
         quad_data, sizeof(quad_data));
  memcpy(solver->quad_adjacency + 4 * solver->max_quad_id,
         adj_data, sizeof(adj_data));

  create_shellable_delta(solver, delta, 0, 0, NULL, 3, quads);
}

void flip_3_cycle(solver *solver, quad_edge edge_a, shellable_delta *delta) {
  quad_edge other_a = *other_side(solver, edge_a);

  quad_edge edge_b = next_edge(other_a);
  quad_edge other_b = *other_side(solver, edge_b);

  quad_edge edge_c = next_edge(other_b);

  quad_edge edges[] = {edge_a, edge_b, edge_c};

  quad_index quads[6] = {
    quad_edge_quad(edge_a),
    quad_edge_quad(edge_b),
    quad_edge_quad(edge_c),

    solver->max_quad_id + 0,
    solver->max_quad_id + 1,
    solver->max_quad_id + 2,
  };

  vertex_index hex[8] = {
    start_v(solver, edge_a),
    start_v(solver, next_edge(edge_a)),
    start_v(solver, opposite_edge(edge_a)),
    start_v(solver, prev_edge(edge_a)),

    end_v(solver, edge_b),
    start_v(solver, opposite_edge(edge_b)),
    solver->max_vertex_id,
    start_v(solver, opposite_edge(edge_c)),
  };

  vertex_index quad_data[3][4];
  quad_edge adj_data[3][4];

  vertex_index finished[1] = {start_v(solver, edge_a)};

  for (size_t i = 0; i < 3; i++) {
    quad_data[i][0] = solver->max_vertex_id;
    quad_data[i][1] = start_v(solver, opposite_edge(edges[i]));
    quad_data[i][2] = start_v(solver, prev_edge(edges[i]));
    quad_data[i][3] = start_v(solver, opposite_edge(edges[(i + 2) % 3]));

    adj_data[i][0] = make_quad_edge(quads[3 + (i + 1) % 3], 3);
    adj_data[i][1] = *other_side(solver, opposite_edge(edges[i]));
    adj_data[i][2] = *other_side(solver, next_edge(edges[(i + 2) % 3]));
    adj_data[i][3] = make_quad_edge(quads[3 + (i + 2) % 3], 0);
  }

  memcpy(solver->hexes + 8 * solver->hex_id, hex, sizeof(hex));
  memcpy(solver->quads + 4 * solver->max_quad_id,
         quad_data, sizeof(quad_data));
  memcpy(solver->quad_adjacency + 4 * solver->max_quad_id,
         adj_data, sizeof(adj_data));

  create_shellable_delta(solver, delta, 1, 1, finished, 3, quads);
}

void shellable_dfs(solver *solver) {
  solver->hex_id++;
  solver->bound_mask = 0;

  if (solver->hex_id == solver->num_hexes) {
    solver->hex_id--;
    return;
  }

  size_t num_children = 0;

  shellable_delta delta;
  no_goods_delta symmetry_delta;

  solver_start_branch(solver, &symmetry_delta);
  uint8_t trail[solver_incremental_symmetry_trail_size(solver)];

  incremental_symmetry_data_trail(solver, trail);

  quad_index boundary[6];
  if (allowed_finish_mesh(solver, boundary)) {
    finish_mesh(solver, boundary);

    solver->bound_mask = ~(uint32_t)0;
    if (!incremental_domination_test(solver)) {
      solver->hex_id++;
      solver->cb(solver, solver->hex_id, solver->user_data);
      solver->hex_id--;

      num_children++;
    }
    incremental_symmetry_data_restore(solver, trail);
    solver->bound_mask = 0;

    undo_finish_mesh(solver, boundary);
  }

  uint32_t num_quads = 0;
  quad_index quads[solver->max_quad_id];

  uint8_t v_min = 0, h_min = 0;
  if ((uint32_t)(solver->hex_id + 1) != solver->num_hexes &&
      solver->num_vertices - solver->max_vertex_id >= 0) {
    for (uint32_t i = solver->quad_queue_offset; i < solver->quad_queue_size; i++) {
      quad_index q = solver->quad_queue[i];
      if (solver->num_occurrences[q] == 1) {
        quads[num_quads++] = q;
      }
    }

    if (solver->num_vertices - solver->max_vertex_id >= 4)
      component_lower_bound(solver, quads, num_quads, &v_min, &h_min);
  }

  if ((uint32_t)(solver->hex_id + 1) == solver->num_hexes ||
      v_min > solver->num_vertices - solver->max_vertex_id ||
      h_min > solver->num_hexes - solver->hex_id ||
      (solver->num_vertices - solver->max_vertex_id >= 4 &&
       look_up_upper_bound(solver, quads, num_quads))) {
    solver_undo_no_goods(solver, &symmetry_delta);
    solver->hex_id--;
    return;
  }

  check_known_solutions(solver, quads, num_quads);

  for (uint32_t i = solver->quad_queue_offset;
       i < solver->quad_queue_size; i++) {
    quad_index q = solver->quad_queue[i];
    if (solver->num_occurrences[q] != 0) {
      if (allowed_flip_5_to_1(solver, q)) {
        flip_5_to_1(solver, q, &delta);
        num_children += recurse_shellable_dfs(solver, &delta, &symmetry_delta, trail);
      }
    }
  }

  for (uint32_t i = solver->quad_queue_offset;
       i < solver->quad_queue_size; i++) {
    quad_index q = solver->quad_queue[i];
    if (solver->num_occurrences[q] != 0) {
      for (quad_edge_index j = 0; j < 4; j++) {
        quad_edge edge = make_quad_edge(q, j);

        if (allowed_flip_4_to_2(solver, edge)) {
          flip_4_to_2(solver, edge, &delta);
          num_children += recurse_shellable_dfs(solver, &delta, &symmetry_delta, trail);
        }

        if (allowed_flip_3_cycle(solver, edge)) {
          flip_3_cycle(solver, edge, &delta);
          num_children += recurse_shellable_dfs(solver, &delta, &symmetry_delta, trail);
        }

        if (allowed_flip_3_line(solver, edge)) {
          flip_3_line(solver, edge, &delta);
          num_children += recurse_shellable_dfs(solver, &delta, &symmetry_delta, trail);
        }

        if (allowed_flip_2_to_4(solver, edge)) {
          flip_2_to_4(solver, edge, &delta);
          num_children += recurse_shellable_dfs(solver, &delta, &symmetry_delta, trail);
        }
      }
    }
  }

  for (uint32_t i = solver->quad_queue_offset;
       i < solver->quad_queue_size; i++) {
    quad_index q = solver->quad_queue[i];
    if (solver->num_occurrences[q] != 0) {
      if (allowed_flip_1_to_5(solver, q)) {
        flip_1_to_5(solver, q, &delta);
        num_children += recurse_shellable_dfs(solver, &delta, &symmetry_delta, trail);
      }
    }
  }

  solver_undo_no_goods(solver, &symmetry_delta);

  if (num_children > 1) solver->num_branches++;

  solver->hex_id--;
}
